#include <src/pcps.h>
#include <src/input.h>
#include <src/correlator.h>
#include <src/wavefunction.h>
#include <src/operator.h>
#include <src/cluster.h>
#include <src/expectation_value.h>
#include <src/rdm.h>
#include <src/blas_lapack.h>

//-------------------------------------------------------------------------------
// pcps::type_lookup -- a function that returns the name of the type whose info
//                      is given
//-------------------------------------------------------------------------------

std::string pcps::type_lookup(const std::type_info & info) {

  std::string retval;

  if      (info == typeid(pcps::OpNN<double>))
    retval = std::string("pcps::OpNN<double>");
  else if (info == typeid(pcps::OpNN<std::complex<double> >))
    retval = std::string("pcps::OpNN<std::complex<double> >");
  else if (info == typeid(pcps::OpRL<double>))
    retval = std::string("pcps::OpRL<double>");
  else if (info == typeid(pcps::OpRL<std::complex<double> >))
    retval = std::string("pcps::OpRL<std::complex<double> >");
  else if (info == typeid(pcps::OpRLRL<double>))
    retval = std::string("pcps::OpRLRL<double>");
  else if (info == typeid(pcps::OpRLRL<std::complex<double> >))
    retval = std::string("pcps::OpRLRL<std::complex<double> >");
  else if (info == typeid(pcps::OpIdentity<double>))
    retval = std::string("pcps::OpIdentity<double>");
  else if (info == typeid(pcps::OpIdentity<std::complex<double> >))
    retval = std::string("pcps::OpIdentity<std::complex<double> >");
  else if (info == typeid(pcps::OpN<double>))
    retval = std::string("pcps::OpN<double>");
  else if (info == typeid(pcps::OpN<std::complex<double> >))
    retval = std::string("pcps::OpN<std::complex<double> >");
  else if (info == typeid(pcps::Part_Fermion))
    retval = std::string("pcps::Part_Fermion");
  else if (info == typeid(pcps::Part_SpinlessFermion))
    retval = std::string("pcps::Part_SpinlessFermion");
  else if (info == typeid(pcps::Part_Spin))
    retval = std::string("pcps::Part_Spin");
  else if (info == typeid(pcps::Ref_RHF))
    retval = std::string("pcps::Ref_RHF");
  else if (info == typeid(pcps::Ref_Uniform))
    retval = std::string("pcps::Ref_Uniform");
  else if (info == typeid(double))
    retval = std::string("double");
  else if (info == typeid(std::complex<double>))
    retval = std::string("std::complex<double>");
  else
    retval = std::string(info.name());

  return retval;

}

//-------------------------------------------------------------------------------
// pcps::recv_expectation_value -- receives an expectation value from a slave
//                                 process and adds it to the appropriate
//                                 element of the ouput array
//-------------------------------------------------------------------------------

template <class SCALAR>
void pcps::recv_expectation_value(int & slave,
                                  int & jobs_in_progress,
                                  SCALAR * const output_array)
{
  // get MPI info
  const MPI::Comm & comm = MPI::COMM_WORLD;
  const int nproc = comm.Get_size();
  const int myrank = comm.Get_rank();

  // prepare a status variable
  MPI::Status status;

  // receive the array position
  int array_pos;
  comm.Recv((void *)(&array_pos), 1, MPI::INT, MPI::ANY_SOURCE, 0, status);

  // determine the slave node that sent the position
  slave = status.Get_source();

  // receive the expectation value
  SCALAR value;
  comm.Recv((void *)(&value), sizeof(SCALAR), MPI::CHAR, slave, 0, status);

  // add the value to the array element
  output_array[array_pos] += value;

  // record the job as finished
  jobs_in_progress--;
}

//-------------------------------------------------------------------------------
// pcps::collect_expectation_values -- collects the expectation values from
//                                     any slave processes that are still working
//                                     and tells all the slaves that there is no
//                                     more work to do
//-------------------------------------------------------------------------------

template <class SCALAR>
void pcps::collect_expectation_values(const std::vector<int> & slaves,
                                      int & jobs_in_progress,
                                      SCALAR * const output_array)
{
  // get MPI info
  const MPI::Comm & comm = MPI::COMM_WORLD;
  const int nproc = comm.Get_size();
  const int myrank = comm.Get_rank();

  // collect expectation values for unfinished jobs
  int slave;
  while (jobs_in_progress > 0)
    pcps::recv_expectation_value(slave, jobs_in_progress, output_array);

  // tell the slaves there is no more work to do
  int signal = -1;
  for (std::vector<int>::const_iterator it = slaves.begin(); it != slaves.end(); it++)
    comm.Send((void *)(&signal), 1, MPI::INT, *it, 0);
}

// instantiate the template
template void pcps::collect_expectation_values(const std::vector<int> &, int &, double               * const);
template void pcps::collect_expectation_values(const std::vector<int> &, int &, std::complex<double> * const);

//-------------------------------------------------------------------------------------
// pcps::expectation_value_detail -- computes the 'projective' expectation value of an
//                                   operator by looping over the configurations of its
//                                   cluster.
//                                   The function is templated on the scalar type
//                                   (SCALAR), the operator type (OP), the particle
//                                   type (PART), and the reference type (REF).
//-------------------------------------------------------------------------------------

template <class SCALAR, class OP, class PART, class REF>
SCALAR pcps::expectation_value_detail(const pcps::Input & userinp,
                                      const pcps::Wavefunction<SCALAR> & wfn,
                                      const OP & op_input,
                                      const int perturbation_level,
                                      pcps::Cluster & local_cluster,
                                      const typename pcps::Correlator<SCALAR>::Set & cluster_corrs,
                                      const typename pcps::Correlator<SCALAR>::Set & pert_cluster_corrs) {

  // get the number of configurations on each lattice site
  const int configs_per_site = userinp.configs_per_site();

  // create a deep copy of the operator to put its data on the top of the memory heap
  const OP op(op_input);

  // organize operator information for fast retrieval during the inner loop
  const int n_op_sites = op.sites().size();
  const int * op_sites = ( n_op_sites > 0 ? &op.sites()[0] : NULL );
  std::vector<int> op_site_ket_occ_vec(n_op_sites);
  int * const op_site_ket_occs = ( n_op_sites > 0 ? &op_site_ket_occ_vec[0] : NULL );

  // organize correlator information for fast retrieval during the inner loop
  const int n_correlators = cluster_corrs.size();
  std::vector<const SCALAR *> corr_data(n_correlators);
  std::vector<const int *> corr_sites(n_correlators);
  std::vector<int> corr_n_sites(n_correlators);
  {
    int i = 0;
    for (typename pcps::Correlator<SCALAR>::Set::const_iterator corr = cluster_corrs.begin(); corr != cluster_corrs.end(); corr++, i++) {
      corr_data[i] = corr->data_ptr();
      corr_n_sites[i] = corr->site_vec().size();
      corr_sites[i] = &corr->site_vec()[0];
    }
  }

  // organize perturbative correlator information for fast retrieval during the inner loop
  const int n_pert_correlators = pert_cluster_corrs.size();
  std::vector<const SCALAR *> pert_corr_data(n_pert_correlators);
  std::vector<const int *> pert_corr_sites(n_pert_correlators);
  std::vector<int> pert_corr_n_sites(n_pert_correlators);
  {
    int i = 0;
    for (typename pcps::Correlator<SCALAR>::Set::const_iterator corr = pert_cluster_corrs.begin(); corr != pert_cluster_corrs.end(); corr++, i++) {
      pert_corr_data[i] = corr->data_ptr();
      pert_corr_n_sites[i] = corr->site_vec().size();
      pert_corr_sites[i] = &corr->site_vec()[0];
    }
  }

  // prepare the function object for computing action values
  pcps::action_value_object<SCALAR, OP, PART> avo(configs_per_site);

  // organize density matrix information for fast retrieval during the inner loop
  pcps::Compact_RDM<SCALAR, OP, PART, REF> rdm(userinp, wfn, op, local_cluster);

  // initialize the expectation value
  SCALAR retval = pcps::zero<SCALAR>();

  // INNER LOOP -- compute the contribution to the expectation value from each configuration of the cluster
  for ( ; !local_cluster.finished(); local_cluster++) {

    SCALAR value = pcps::unity<SCALAR>();

    // calculate the numerator and denominator contributions for each correlator
    for (int i = 0; i < n_correlators; i++) {
      int cmpd_ket_index = 0;
      int cmpd_bra_index = 0;
      for (int j = 0; j < corr_n_sites[i]; j++) {
        cmpd_ket_index = cmpd_ket_index * configs_per_site + local_cluster.ket_occ(corr_sites[i][j]);
        cmpd_bra_index = cmpd_bra_index * configs_per_site + local_cluster.bra_occ(corr_sites[i][j]);
      }
      const SCALAR ket_element = corr_data[i][cmpd_ket_index];
      const SCALAR bra_element = corr_data[i][cmpd_bra_index];
      if ( std::abs(bra_element) > 1.0e-5 )
        value *= ket_element / bra_element;
      else
        value = pcps::zero<SCALAR>();
    }

    // calculate the perturbative correlator contributions
    if (n_pert_correlators == 1) {
      int cmpd_ket_index = 0;
      int cmpd_bra_index = 0;
      for (int j = 0; j < pert_corr_n_sites[0]; j++) {
        cmpd_ket_index = cmpd_ket_index * configs_per_site + local_cluster.ket_occ(pert_corr_sites[0][j]);
        cmpd_bra_index = cmpd_bra_index * configs_per_site + local_cluster.bra_occ(pert_corr_sites[0][j]);
      }
      const SCALAR ket_element = pert_corr_data[0][cmpd_ket_index];
      const SCALAR bra_element = pert_corr_data[0][cmpd_bra_index];
      if (perturbation_level >= 2)
        value *= ket_element - bra_element * ( pcps::unity<SCALAR>() - bra_element + ket_element );
      else
        value *= ket_element - bra_element;
    } else if (n_pert_correlators == 2) {
      int cmpd_ket_index = 0;
      int cmpd_bra_index = 0;
      for (int j = 0; j < pert_corr_n_sites[0]; j++) {
        cmpd_ket_index = cmpd_ket_index * configs_per_site + local_cluster.ket_occ(pert_corr_sites[0][j]);
        cmpd_bra_index = cmpd_bra_index * configs_per_site + local_cluster.bra_occ(pert_corr_sites[0][j]);
      }
      const SCALAR ket_element_0 = pert_corr_data[0][cmpd_ket_index];
      const SCALAR bra_element_0 = pert_corr_data[0][cmpd_bra_index];
      cmpd_ket_index = 0;
      cmpd_bra_index = 0;
      for (int j = 0; j < pert_corr_n_sites[1]; j++) {
        cmpd_ket_index = cmpd_ket_index * configs_per_site + local_cluster.ket_occ(pert_corr_sites[1][j]);
        cmpd_bra_index = cmpd_bra_index * configs_per_site + local_cluster.bra_occ(pert_corr_sites[1][j]);
      }
      const SCALAR ket_element_1 = pert_corr_data[1][cmpd_ket_index];
      const SCALAR bra_element_1 = pert_corr_data[1][cmpd_bra_index];
      value *= ( ket_element_0 - bra_element_0 ) * ( ket_element_1 - bra_element_1 );
    }

    // calculate the contribution of the operator's action
    for (int i = 0; i < n_op_sites; i++)
      op_site_ket_occs[i] = local_cluster.ket_occ(op_sites[i]);
    value *= avo.get_value(op_site_ket_occs, configs_per_site);

    // calculate the contribution from the density matrix element
    value *= rdm.element(local_cluster);

    // add the total contribution to the expectation value
    retval += value;

  }
  // END OF INNER LOOP

  // multiply the expectation value by the operator's numeric coefficient
  retval *= op.coeff();

  // return the expectation value
  return retval;

}

//-------------------------------------------------------------------------------------
// pcps::expectation_value_choose_ref -- Determines which wavefunction reference is
//                                       being used and calls the appropriate template.
//
//                                       This function is templated on the scalar type
//                                       (SCALAR), the operator type (OP), and the
//                                       particle type (PART).
//-------------------------------------------------------------------------------------

template <class SCALAR, class OP, class PART>
SCALAR pcps::expectation_value_choose_ref(const pcps::Input & userinp,
                                          const pcps::Wavefunction<SCALAR> & wfn,
                                          const OP & op,
                                          const int perturbation_level,
                                          pcps::Cluster & local_cluster,
                                          const typename pcps::Correlator<SCALAR>::Set & cluster_corrs,
                                          const typename pcps::Correlator<SCALAR>::Set & pert_cluster_corrs) {

  // function for uniform reference
  if ( userinp.ref_type() == pcps::keyword::Ref_Uniform )
    return pcps::expectation_value_detail<SCALAR, OP, PART, pcps::Ref_Uniform>
                                   (userinp, wfn, op, perturbation_level, local_cluster, cluster_corrs, pert_cluster_corrs);

  // function for rhf reference
  else if ( userinp.ref_type() == pcps::keyword::Ref_RHF )
    return pcps::expectation_value_detail<SCALAR, OP, PART, pcps::Ref_RHF>
                                   (userinp, wfn, op, perturbation_level, local_cluster, cluster_corrs, pert_cluster_corrs);

  // unknown reference type
  throw pcps::Exception( (boost::format("expectation_value encountered unexpected reference type \"%s\"") % userinp.ref_type()).str() );
  return pcps::zero<SCALAR>();

}

//-------------------------------------------------------------------------------------
// pcps::expectation_value_choose_part -- Determines which particle type is being used
//                                        and calls the appropriate template.
//
//                                        This function is templated on the scalar type
//                                        (SCALAR) and the operator type (OP).
//-------------------------------------------------------------------------------------

template <class SCALAR, class OP>
SCALAR pcps::expectation_value_choose_part(const pcps::Input & userinp,
                                           const pcps::Wavefunction<SCALAR> & wfn,
                                           const OP & op,
                                           const int perturbation_level,
                                           pcps::Cluster & local_cluster,
                                           const typename pcps::Correlator<SCALAR>::Set & cluster_corrs,
                                           const typename pcps::Correlator<SCALAR>::Set & pert_cluster_corrs) {

  // function for fermions
  if ( userinp.particle_type() == pcps::keyword::Part_Fermion )
    return pcps::expectation_value_choose_ref<SCALAR, OP, pcps::Part_Fermion>
                                             (userinp, wfn, op, perturbation_level, local_cluster, cluster_corrs, pert_cluster_corrs);

  // function for spinless fermions
  else if ( userinp.particle_type() == pcps::keyword::Part_Spinless_Fermion )
    return pcps::expectation_value_choose_ref<SCALAR, OP, pcps::Part_SpinlessFermion>
                                             (userinp, wfn, op, perturbation_level, local_cluster, cluster_corrs, pert_cluster_corrs);

  // function for spins
  else if ( userinp.particle_type() == pcps::keyword::Part_Spin )
    return pcps::expectation_value_choose_ref<SCALAR, OP, pcps::Part_Spin>
                                             (userinp, wfn, op, perturbation_level, local_cluster, cluster_corrs, pert_cluster_corrs);

  // unknown particle type
  throw pcps::Exception( (boost::format("expectation_value encountered unexpected particle type \"%s\"") % userinp.particle_type()).str() );
  return pcps::zero<SCALAR>();

}

//-------------------------------------------------------------------------------------
// pcps::expectation_value_choose_op -- Determines which operator type is being used
//                                      and calls the appropriate template.
//
//                                      This function is templated on the scalar
//                                      type (SCALAR).
//-------------------------------------------------------------------------------------

template <class SCALAR>
SCALAR pcps::expectation_value_choose_op(const pcps::Input & userinp,
                                         const pcps::Wavefunction<SCALAR> & wfn,
                                         const typename pcps::OpBase<SCALAR>::const_ptr & op,
                                         const int perturbation_level,
                                         pcps::Cluster & local_cluster,
                                         const typename pcps::Correlator<SCALAR>::Set & cluster_corrs,
                                         const typename pcps::Correlator<SCALAR>::Set & pert_cluster_corrs) {

  // function for raising/lowering operator
  if ( typeid(*op) == typeid(pcps::OpRL<SCALAR>) ) {
    boost::shared_ptr<const pcps::OpRL<SCALAR> > derived_op_ptr = boost::dynamic_pointer_cast<const pcps::OpRL<SCALAR>, const pcps::OpBase<SCALAR> >(op);
    return pcps::expectation_value_choose_part<SCALAR, pcps::OpRL<SCALAR> >
                                              (userinp, wfn, *derived_op_ptr, perturbation_level, local_cluster, cluster_corrs, pert_cluster_corrs);
  }

  // function for double number operator
  else if ( typeid(*op) == typeid(pcps::OpNN<SCALAR>) ) {
    boost::shared_ptr<const pcps::OpNN<SCALAR> > derived_op_ptr = boost::dynamic_pointer_cast<const pcps::OpNN<SCALAR>, const pcps::OpBase<SCALAR> >(op);
    return pcps::expectation_value_choose_part<SCALAR, pcps::OpNN<SCALAR> >
                                              (userinp, wfn, *derived_op_ptr, perturbation_level, local_cluster, cluster_corrs, pert_cluster_corrs);
  }

  // function for single number operator
  else if ( typeid(*op) == typeid(pcps::OpN<SCALAR>) ) {
    boost::shared_ptr<const pcps::OpN<SCALAR> > derived_op_ptr = boost::dynamic_pointer_cast<const pcps::OpN<SCALAR>, const pcps::OpBase<SCALAR> >(op);
    return pcps::expectation_value_choose_part<SCALAR, pcps::OpN<SCALAR> >
                                              (userinp, wfn, *derived_op_ptr, perturbation_level, local_cluster, cluster_corrs, pert_cluster_corrs);
  }

  // function for identity operator
  else if ( typeid(*op) == typeid(pcps::OpIdentity<SCALAR>) ) {
    boost::shared_ptr<const pcps::OpIdentity<SCALAR> > derived_op_ptr = boost::dynamic_pointer_cast<const pcps::OpIdentity<SCALAR>, const pcps::OpBase<SCALAR> >(op);
    return pcps::expectation_value_choose_part<SCALAR, pcps::OpIdentity<SCALAR> >
                                              (userinp, wfn, *derived_op_ptr, perturbation_level, local_cluster, cluster_corrs, pert_cluster_corrs);
  }

  // unknown operator type
  throw pcps::Exception( ( boost::format("expectation_value encountered unexpected operator of type \"%s\"") % pcps::type_lookup(typeid(*op)) ).str() );
  return pcps::zero<SCALAR>();

}

//-------------------------------------------------------------------------------
// pcps::expectation_value_master -- assigns slaves to compute an expectation
//                                   value, or computes it directly if there are
//                                   no slaves
//-------------------------------------------------------------------------------

template <class SCALAR>
void pcps::expectation_value_master(const std::vector<int> & slaves,
                                    int & slaves_used,
                                    int & jobs_in_progress,
                                    SCALAR * const output_array,
                                    const int array_pos,
                                    const SCALAR coeff,
                                    const pcps::Input & userinp,
                                    const pcps::Wavefunction<SCALAR> & wfn,
                                    const typename pcps::OpBase<SCALAR>::const_ptr & op,
                                    const int perturbation_level,
                                    pcps::Cluster & local_cluster,
                                    const std::vector<int> & corr_pos,
                                    const std::vector<int> & pert_corr_pos) {

  // get MPI info
  const MPI::Comm & comm = MPI::COMM_WORLD;
  const int nproc = comm.Get_size();
  const int myrank = comm.Get_rank();

  // if there are no slaves or we are in site parallel mode, compute the expectation directly
  if (slaves.size() == 0 || userinp.parallel_type() == pcps::keyword::Parallel_Site) {

    // get the set of non-perturbative correlators
    typename pcps::Correlator<SCALAR>::Set corrs;
    for (std::vector<int>::const_iterator pos = corr_pos.begin(); pos != corr_pos.end(); pos++)
      corrs.insert(wfn.correlators().at(*pos));

    // get the set of perturbartive correlators
    typename pcps::Correlator<SCALAR>::Set pert_corrs;
    for (std::vector<int>::const_iterator pos = pert_corr_pos.begin(); pos != pert_corr_pos.end(); pos++)
      pert_corrs.insert(wfn.pert_corrs().at(*pos));

    // compute the expectation value
    output_array[array_pos] += coeff * expectation_value_choose_op<SCALAR>(userinp, wfn, op, perturbation_level, local_cluster, corrs, pert_corrs);
    return;

  }

  // split the cluster into multiple clusters such that each of the new clusters has fewer than max_configs configurations
  const int max_configs = 5000000;
  std::list<pcps::Cluster> subclusters;
  local_cluster.split(max_configs, subclusters);

  // assign work to slaves
  for (std::list<pcps::Cluster>::const_iterator clust = subclusters.begin(); clust != subclusters.end(); clust++) {

    int slave;

    // if all the slaves are busy, wait for one to finish
    if (slaves_used == slaves.size())
      pcps::recv_expectation_value(slave, jobs_in_progress, output_array);

    // if not all the slaves are in use, find one that is not working
    else if (slaves_used < slaves.size()) {
      slave = slaves.at(slaves_used);
      slaves_used++;
    }

    // raise an error if we have counted the slaves incorrectly
    else
      throw pcps::Exception("slaves_used is greater than the number of slaves");

    // send a signal telling the slave it has work to do
    int signal = 1;
    comm.Send((void *)(&signal), 1, MPI::INT, slave, 0);

    // send the operator to the slave
    pcps::send_operator<SCALAR>(comm, op, slave, 0);

    // send the perturbation level to the slave
    comm.Send((const void *)(&perturbation_level), 1, MPI::INT, slave, 0);

    // send the cluster to the slave
    clust->send(comm, slave, 0);

    // send the non-perturbative correlator positions to the slave
    pcps::send_vec(comm, corr_pos, slave, 0);

    // send the perturbative correlator positions to the slave
    pcps::send_vec(comm, pert_corr_pos, slave, 0);

    // send the array position to the slave
    comm.Send((const void *)(&array_pos), 1, MPI::INT, slave, 0);

    // send the numeric coefficient to the slave
    comm.Send((const void *)(&coeff), sizeof(SCALAR), MPI::CHAR, slave, 0);

    // count the job
    jobs_in_progress++;

  }

}

//-------------------------------------------------------------------------------
// pcps::expectation_value_slave -- places a slave process in a state in which it
//                                  receives instructions to compute expectation
//                                  values from the master process
//-------------------------------------------------------------------------------

template <class SCALAR>
void pcps::expectation_value_slave(const pcps::Input & userinp,
                                   const pcps::Wavefunction<SCALAR> & wfn) {

  // get MPI info
  const MPI::Comm & comm = MPI::COMM_WORLD;
  const int nproc = comm.Get_size();
  const int myrank = comm.Get_rank();

  // do work as instructed by the master process
  while (true) {

    // receive the signal telling us whether there is more work to do
    int signal;
    comm.Recv((void *)(&signal), 1, MPI::INT, 0, 0);
    if (signal < 0) break;

    int array_pos;
    SCALAR value;
    {

      // receive the operator
      typename pcps::OpBase<SCALAR>::const_ptr op;
      pcps::recv_operator<SCALAR>(comm, op, 0, 0);

      // receive the perturbation level
      int perturbation_level;
      comm.Recv((void *)(&perturbation_level), 1, MPI::INT, 0, 0);

      // receive the cluster information
      pcps::Cluster local_cluster;
      local_cluster.recv(comm, 0, 0);

      // receive the non-perturbative correlator positions
      std::vector<int> corr_pos;
      pcps::recv_vec(comm, corr_pos, 0, 0);

      // receive the perturbative correlator positions
      std::vector<int> pert_corr_pos;
      pcps::recv_vec(comm, pert_corr_pos, 0, 0);

      // receive the array position
      comm.Recv((void *)(&array_pos), 1, MPI::INT, 0, 0);

      // receive the numeric coefficient
      SCALAR coeff;
      comm.Recv((void *)(&coeff), sizeof(SCALAR), MPI::CHAR, 0, 0);

      // get the set of non-perturbative correlators
      typename pcps::Correlator<SCALAR>::Set corrs;
      for (std::vector<int>::const_iterator pos = corr_pos.begin(); pos != corr_pos.end(); pos++)
        corrs.insert(wfn.correlators().at(*pos));

      // get the set of perturbartive correlators
      typename pcps::Correlator<SCALAR>::Set pert_corrs;
      for (std::vector<int>::const_iterator pos = pert_corr_pos.begin(); pos != pert_corr_pos.end(); pos++)
        pert_corrs.insert(wfn.pert_corrs().at(*pos));

      // compute the expectation value
      value = coeff * expectation_value_choose_op<SCALAR>(userinp, wfn, op, perturbation_level, local_cluster, corrs, pert_corrs);

    }

    // send the array position back to the master process
    comm.Send((void *)(&array_pos), 1, MPI::INT, 0, 0);

    // send the value back to the master process
    comm.Send((const void *)(&value), sizeof(SCALAR), MPI::CHAR, 0, 0);

  }

}

// instantiate the template
template void pcps::expectation_value_slave(const pcps::Input &, const pcps::Wavefunction< double               > &);
template void pcps::expectation_value_slave(const pcps::Input &, const pcps::Wavefunction< std::complex<double> > &);

//-------------------------------------------------------------------------------
// pcps::expectation_value -- Computes the projective cps approximation to the
//                            operator's expectation value.
//                            The operator is stored in op.
//                            The expectation value is computed up to order
//                            perturbation_level in the perturbative correlators.
//
//                            This particular function constructs the cluster
//                            around the operator and determines the relavant
//                            correlators.
//                            The expectation_value_master function is then
//                            called, which deals with computing the result in
//                            parallel if farm parallel mode has been specified.
//                            Eventually, after a number of other function calls
//                            to determine the operator, particle, and reference
//                            types, the expectation_value_detail function is
//                            called, which computes the expectation value.
//-------------------------------------------------------------------------------

template<class SCALAR>
void pcps::expectation_value(const std::vector<int> & slaves,
                             int & slaves_used,
                             int & jobs_in_progress,
                             SCALAR * const output_array,
                             const int array_pos,
                             const SCALAR coeff,
                             const pcps::Input & userinp,
                             const pcps::Wavefunction<SCALAR> & wfn,
                             const typename pcps::OpBase<SCALAR>::const_ptr & op,
                             const int perturbation_level,
                             const pcps::Cluster & fixed_bra_occs,
                             const pcps::Cluster & fixed_ket_occs,
                             const typename pcps::Correlator<SCALAR>::Set & excluded) {

  // check that the perturbation level is implemented
  if (perturbation_level < 0 || perturbation_level > 2)
    throw pcps::Exception("perturbation_level must be 0, 1, or 2");

  // if the cluster is empty and the operator is the identity, the expectation value is trivial
  if ( typeid(*op) == typeid(pcps::OpIdentity<SCALAR>) && fixed_bra_occs.size() == 0 && fixed_ket_occs.size() == 0 ) {
    output_array[array_pos] += coeff * op->coeff();
    return;
  }

  // get the non-perturbative correlators involved in the cluster
  std::vector<int> corr_pos;
  {
    typename pcps::Correlator<SCALAR>::Set cluster_corrs;
    for (int i = 0; i < op->sites().size(); i++)
      if (op->offsets().at(i) != 0) {
        //const std::vector<pcps::Correlator> & touching = wfn.correlators_with_site(op->sites().at(i));
        //for (std::vector<pcps::Correlator>::const_iterator corr = touching.begin(); corr != touching.end(); corr++)
        //  if (excluded.count(*corr) == 0 && cluster_corrs.count(*corr) == 0)
        for (typename std::vector< pcps::Correlator<SCALAR> >::const_iterator corr = wfn.correlators().begin(); corr != wfn.correlators().end(); corr++)
          if (    std::count(corr->site_vec().begin(), corr->site_vec().end(), op->sites().at(i)) > 0
               && excluded.count(*corr) == 0
               && cluster_corrs.count(*corr) == 0 )
            cluster_corrs.insert(*corr); // shallow copy
      }
    for (int i = 0; i < wfn.correlators().size(); i++)
      if ( cluster_corrs.count(wfn.correlators().at(i)) > 0 )
        corr_pos.push_back(i);
  }

  // create a vector to hold the positions of perturbative correlators that are currently in use
  std::vector<int> pert_corr_pos;

  // compute the non-pertrubative expectation value
  if (perturbation_level >= 0) {

    // initialize the cluster of sites
    pcps::Cluster local_cluster(userinp, wfn, op, fixed_bra_occs, fixed_ket_occs);

    // compute the cluster's non-perturbative expectation value
    if (local_cluster.nconfigs() > 0)
      pcps::expectation_value_master<SCALAR>(slaves, slaves_used, jobs_in_progress, output_array, array_pos, coeff,
                                             userinp, wfn, op, perturbation_level, local_cluster, corr_pos, pert_corr_pos);

  }

  // declare a vector to hold the positions of the perturbative correlators affected by the operator
  std::vector<int> affected_pert_corr_pos;

  // compute the 1st order pertrubative correction
  if (perturbation_level >= 1) {

    // populate the set of perturbative correlators affected by the operator
    {
      typename pcps::Correlator<SCALAR>::Set affected_pert_corrs;
      for (int i = 0; i < op->sites().size(); i++)
        if (op->offsets().at(i) != 0)
          for (typename std::vector< pcps::Correlator<SCALAR> >::const_iterator corr = wfn.pert_corrs().begin(); corr != wfn.pert_corrs().end(); corr++)
            if (    std::count(corr->site_vec().begin(), corr->site_vec().end(), op->sites().at(i) ) > 0
                 && excluded.count(*corr) == 0
                 && affected_pert_corrs.count(*corr) == 0 )
              affected_pert_corrs.insert(*corr); // shallow copy
      for (int i = 0; i < wfn.pert_corrs().size(); i++)
        if ( affected_pert_corrs.count(wfn.pert_corrs().at(i)) > 0 )
          affected_pert_corr_pos.push_back(i);
    }

    // compute the correction for each perturbative correlator affected by the operator
    for (std::vector<int>::const_iterator pos = affected_pert_corr_pos.begin(); pos != affected_pert_corr_pos.end(); pos++) {

      // place the perturbative correlator in the set of "in use" perturbative correlators
      pert_corr_pos.clear();
      pert_corr_pos.push_back(*pos);

      // initialize the cluster of sites
      pcps::Cluster local_cluster(userinp, wfn, op, fixed_bra_occs, fixed_ket_occs, pert_corr_pos);

      // compute the cluster's perturbative expectation value
      if (local_cluster.nconfigs() > 0)
        pcps::expectation_value_master<SCALAR>(slaves, slaves_used, jobs_in_progress, output_array, array_pos, coeff,
                                               userinp, wfn, op, perturbation_level, local_cluster, corr_pos, pert_corr_pos);

    }

  }

  // compute the 2nd order pertrubative correction
  if (perturbation_level >= 2) {

    // compute the correction for each unique pair of affected perturbative correlators
    for (std::vector<int>::const_iterator pos1 = affected_pert_corr_pos.begin(); pos1 != affected_pert_corr_pos.end(); pos1++) {
      std::vector<int>::const_iterator pos2 = pos1;
      pos2++;
      for ( ; pos2 != affected_pert_corr_pos.end(); pos2++) {

        // place the perturbative correlator in the set of "in use" perturbative correlators
        pert_corr_pos.clear();
        pert_corr_pos.push_back(*pos1);
        pert_corr_pos.push_back(*pos2);

        // initialize the cluster of sites
        pcps::Cluster local_cluster(userinp, wfn, op, fixed_bra_occs, fixed_ket_occs, pert_corr_pos);

        // compute the cluster's perturbative expectation value
        if (local_cluster.nconfigs() > 0)
          pcps::expectation_value_master<SCALAR>(slaves, slaves_used, jobs_in_progress, output_array, array_pos, coeff,
                                                 userinp, wfn, op, perturbation_level, local_cluster, corr_pos, pert_corr_pos);

      }
    }

  }

}

// instantiate the template for doubles
template void pcps::expectation_value(const std::vector<int> & slaves,
                                      int & slaves_used,
                                      int & jobs_in_progress,
                                      double * const output_array,
                                      const int array_pos,
                                      const double coeff,
                                      const pcps::Input & userinp,
                                      const pcps::Wavefunction<double> & wfn,
                                      const pcps::OpBase<double>::const_ptr & op,
                                      const int perturbation_level,
                                      const pcps::Cluster & fixed_bra_occs,
                                      const pcps::Cluster & fixed_ket_occs,
                                      const pcps::Correlator<double>::Set & excluded);

// instantiate the template for complex doubles
template void pcps::expectation_value(const std::vector<int> & slaves,
                                      int & slaves_used,
                                      int & jobs_in_progress,
                                      std::complex<double> * const output_array,
                                      const int array_pos,
                                      const std::complex<double> coeff,
                                      const pcps::Input & userinp,
                                      const pcps::Wavefunction< std::complex<double> > & wfn,
                                      const pcps::OpBase<std::complex<double> >::const_ptr & op,
                                      const int perturbation_level,
                                      const pcps::Cluster & fixed_bra_occs,
                                      const pcps::Cluster & fixed_ket_occs,
                                      const pcps::Correlator< std::complex<double> >::Set & excluded);
